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Abstract 

Computational time reversal imaging can be used to locate the position 
of multiple scatterers in a known background medium. Here, we discuss 
a sparse approximation method for computational time-reversal imaging. 
The method is formulated entirely in the frequency domain, and besides 
imaging it can also be used for denoising, and to determine the magnitude 
of the scattering coefficients in the presence of moderate noise levels. 

PACS: 

02.30.Zz Inverse problems 
43.60.+d Acoustic signal processing 

43.60.Pt Signal processing techniques for acoustic inverse problems 
43.60.Tj Wave front reconstruction, acoustic time-reversal 
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1 Introduction 

Acoustic, elastic or electro-magnetic waves scattered by the inhomogeneities in 
the medium carry significant information, which can be used to obtain true 
images of the investigated domain [1]. Recently, it was shown that scattered 
acoustic waves can be time-reversed and focused onto their original source lo- 
cation through arbitrary media, using a so-called time- reversal mirror [2] . This 
important result shows that one can use computational time reversal imaging 
to identify the location of multiple point scatterers (targets) in a known back- 
ground medium [3]. In this back-propagated signal is computed, rather 
than implemented in the real medium, and its peaks indicate the existence 
of possible scattering targets. The current algorithms for computational time 
reversal imaging are based on the null subspace projection operator, obtained 
through the singular value decomposition of the frequency response matrix [4-9] . 
Here, we discuss a different approach based on a sparse approximation method. 
Besides imaging, this approach can be used for denoising, and to determine the 
magnitude of the scattering coefficients of the targets embedded in homogeneous 
media in the presence of moderate noise levels. The method is formulated en- 
tirely in the frequency domain, in the case where the Born approximation can 
be used [10]. Relevant applications are radar imaging, exploration seismics, 
nondestructive material testing, microwave breast imaging, ultrasound kidney 
stones localization, and other acoustic inverse problems [1-9]. 

2 Frequency response matrix 

We consider a system consisting of an array of N transceivers (i.e. each antenna 
is an emitter and a receiver) located at x n £ R D (n = 1, N), and a collection 
of M distinct scatterers (targets) with scattering coefficients p m , located at 
y m £ R D (m = 1, M) (Fig. 1). Here, D = 1, 2, 3 is the dimensionality of the 
space. Also, we assume that the wave propagation is well approximated in the 
space-frequency domain (x, u) by the inhomogeneous Helmholtz equation: 



where ip(x, uS) is the wave amplitude produced by a localized source s(x, w), ko = 
2ttui/c = 2n/X is the wavenumber of the homogeneous background, with uu the 
frequency, Co the homogeneous background wave speed, and A the wavelength. 
Here, rj(x) is the index of refraction: r](x) = cq/c(x), where c(x) is the wave 
speed at location x. In the background we have tJq{x) = 1, while if{x) = 
1 + a(x), measures the change in the wave speed at the scatterers location. 

The fundamental solutions, or the Green functions, for this problem satisfy 
the following equations: 



[V 2 + k 2 r] 2 (x)] ip{x,uj) = -s(x,oj), 



(1) 



[V 2 + k 2 ] G {x,x') = -S(x-x'), 



(2) 




(3) 
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for the homogeneous and inhomogencous media, respectively. The expression 
of the homogeneous Green function depends on the dimension of the space as 
following: 

G (x,x') = —exp(ik \x-x'\), D = 1, (4) 

G (x,x') = l -H^\k \x-x'\), D = 2, (5) 
exp(ik \x -a;'|) 

G (x,x) = , D = 3, (6) 

where, H^p {.) is the zero order Hankel function of the first kind. Also, the 
fundamental solution G(x, x') for the inhomogeneous medium can be written in 
terms of that for the homogeneous one G (x,x') as: 

G(x, x') = G (x, x) + jfeg J a{z)G {x, z)G{z, x')dz. (7) 

This is an implicit integral equation for G(x,x'). Since the scatterers are as- 
sumed to be pointlike, the regions with a(z) ^ are assumed to be finite, and 
included in compact domains Sl m centered at y m , m = 1, M, which are small 
compared to the wavelength A. Therefore we can write: 

M 

a(z,w) = ^2 Pm(u)S(z - y m ), (8) 

m— 1 

and consequently we obtain the Foldy-Lax equations: 

M 

G(x, x') ~ G {x, x') + ^2 p m {u))G (x, y m )G(y m , x'). (9) 

m— 1 

If the scatterers are sufficiently far apart then we can neglect the multiple scat- 
tering among the scatterers (G(y m ,x') ~ Go(y m ,x')) and we obtain the Born 
approximation of the solution: 

M 

G(x,x') ~ G (a;, a;') + ^ p m (ui)G (x, y m )G {y m , x'). (10) 

m— 1 

If x corresponds to the receiver location Xi, and x' corresponds to the emitter 
location Xj, then we obtain: 

G(xi,Xj) ~ G (xi,Xj) +Hij((jj), (11) 

where 

M 

Hiji^) = ^2 G o(xi,y m )Pm(u)Go(y m ,Xj), i,j = 1, N, (12) 

m—l 

are the elements of the frequency response matrix H(lu) = [Hij(u>)] [4-9]. The 
response matrix H(oS) is obviously a complex and symmetric N x N matrix, 
since the same Green function is used in both the transmission and the reception 
paths. 
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3 Back-propagation and imaging 

An important step in computational time-reversal imaging is to determine the 
frequency response matrix H{lo). This can be done by performing a series of 
N simple experiments, in which a single clement of the array is excited with a 
suitable signal s and we measure the frequency response between this element 
and all the other elements of the array [1-9]. In general, given the Green func- 
tion of the homogeneous media Gq(x, x'), the general solution to the Hclmholtz 
equation is the convolution: 

i/>(x, u) = (Go * s)(x) = J G {x, x')s(x', u)dx' . (13) 

Thus, if the j antenna emits a signal Sj then, using the convolution theorem 
in the Fourier domain, the field produced at the location r is Go(r,Xj)sj. If 
this field is incident on the m-th scatterer, it produces at r the scattered field 
Go(r,y m )p m (u))Go(r,Xj)sj. Thus, the total wave field at location r, due to a 
pulse emitted by a single element at xj and scattered by the M targets can be 
expressed as: 

M 

ip{r,u) = ^2 G (r,y m )p m (w)G (y m ,Xj)sj. (14) 

m— 1 

If this field is measured at the i-th antenna we obtain: 

M 

ip(xi,u) = ^2 G (xi,y m )p m (u>)G (y m ,Xj)sj = Hij(u>)sj. (15) 

m— 1 

Thus, the response matrix can be rewritten as: 

H(u) = r(oj)R(uj)r T (oj), (16) 

where 

R(u>) = diag{pi{uj), ...,p M (w)}, (17) 

and 

r(w) =[ g(yi, w) g(yi,u) ... g{y M ,u)]. (18) 
Here, g(y m ,w) is the Green function vector associated with the m-th scatterer 

g(y m ,w) = [ G {xi,y m ,u) G (x 2 ,y m ,w) ... G (x N ,y m , w) ] T . (19) 

In general, the Green function vector defined as: 

g{r,u) = [ G (xi,r,ui) G (x 2 ,r,uj) ... G a (x N , r, d) ] T , (20) 

expresses the response at each array element due to a single pulse emitted from 
r. 

In the formulation of time-reversal imaging one forms the self-adjoint matrix 

[6]: 

K(w) = H*(u)H(uj) = H(uj)H{lj), (21) 
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where the star denotes the adjoint and the bar denotes the complex conjugate 
(H* = H, since H is symmetric). H is the frequency-domain version of a time- 
reversed response matrix, thus K(w) corresponds to performing a scattering 
experiment, time-reversing the received signals and using them as input for a 
second scattering experiment. Therefore, time- reversal imaging relies on the 
assumption that the Green function can be always calculated. 

As long as the number of antenna elements exceeds the number of scatter- 
ers, M < N, the matrix K(u>) is rank deficient and it has only M non-zero 
eigenvalues, with the corresponding eigenvectors V m (u), m = 1,...,M. When 
the scatterers are well resolved, the columns of the matrix T(ui) are approxi- 
mately orthogonal to each other, and the eigenvectors can be back-propagated 
as g T (r,ui)V m (uj), and consequently the radiated wavefields focus at target lo- 
cations. Thus, each eigenvector can be used to locate a single scatterer. For 
example, let us consider a scenario consisting of N = 100 transceivers, sepa- 
rated by d = A/2, and located at x n = [ n\/2 + a/2 - NX/4 ] T . Also, 
there are two targets, M = 2, with the scattering coefficients p 12 = 1, situ- 
ated at yi = [ 0.65 0.25 ] T a, and respectively tj2 = [ 0.80 0.75 ] T a. Here, 
a = 100A is the side of the imaging area, and the computational image grid is set 
to L x L = 200 x 200 pixels. In Figure 2 we give the first two back-propagated 
eigenvectors and the computed time-reversal image. One can see that since the 
targets are well separated the computed time-reversal image is almost a perfect 
superposition of the two independently back-propagated eigenvectors. 

4 Subspace-based imaging 

The above result does not apply to the case of poorly-resolved targets. In 
this case, the eigenvectors of K(u>) are linear combinations of the target Green 
function vectors g{y m ,uj). Thus, back-propagating one of these eigenvectors 
generates a linear combination of wavefields, each focused on a different target 
location. The subspace-based algorithms, based on the multiple signal classifi- 
cation (MUSIC) method, can be used in this more general situation [7-9]. The 
signal subspace method assumes that the number M of point targets in the 
medium is lower than the number of transceivers TV, and the general idea is to 
localize multiple sources by exploiting the eigenstructure and the rank deficiency 
of the response matrix H(u>). 

The singular value decomposition of the symmetrical matrix H(u) is given 



where U (u>) and V(u)) are the N x N orthogonal matrices corresponding to the 
left and right singular vectors (V(co) = U(u), since H(u) is symmetric). The 
singular value matrix A(w) is diagonal A(o>) = diag{X\(u)), Aat(w)}, and since 
H(u) is rank-deficient, all but the first M singular values vanish: 



by: 



H(u) = U(uj)A(uj)V*(uj), 



(22) 



Xi(u>) ± 0, i = l,...,M, 
A» = 0, j = M + l,...,N. 



(23) 
(24) 



5 



Therefore, the first M columns (singular vectors) of V(u>) span the same sub- 
space a as the columns of r(w), while the remaining N — M columns span the 
null-subspacc v of T(uj). Thus, by partitioning V(u) as: 

V(w) = [ KM KM ] (25) 

where V a {u) has the first M columns and V v (w) has the remaining N — M 
columns, one can write the signal space as a direct sum a®i>, where the essential 
signal-subspace a is orthogonal to the null-subspace v. It follows immediately 
that: 

K»ff(w) = 0, (26) 

and therefore 

TOs(r,«) = 0, (27) 

for any u>. Therefore, the target locations must correspond to the peaks in the 
MUSIC pseudo-spectrum for any w. 

P M usic(r,u) = \\V:(uj)g(r,uj)\\- 2 , (28) 

where g(r,u>) is the free-space Green function vector. Thus, one can form an 
image of the scatterers by plotting, at each point r, the quantity PMUSic{ r , LJ )- 
The resulting plot will have large peaks at the locations of the scatterers. For 
example, let us consider a two dimensional scenario, consisting of N = 100 
linearly distributed transceivers, separated by d = A/2 and located at x n = 
[ n\/2 + a/2 — NX/ 4 ] , where a = 100 A is the side of the imaging area. 
Also, there are M = 5 targets with the scattering coefficients p m = 10, m = 
1, M. The computational image grid is set to L x L = 200 x 200 pixels. We 
consider two cases, one without noise, and the second one with Gaussian noise 
added to the elements of the response matrix Hij(w). The noise level was set 
such that the signal to noise ratio (SNR) is SNR — 2. SNR compares the level 
of a desired signal to the level of background noise. The higher the ratio, the 
less obtrusive the background noise is. SNR measures the power ratio between 
a signal and the background noise: 

SNR = P signal / Pnoise — {^-signal / ^-noise ) : (^) 

where P is average power and A is root mean square (RMS) amplitude. In 
Figure 3 we give the results obtained with the MUSIC algorithm in both cases. 
One can see that the MUSIC pseudo-spectrum provides a better resolution and 
separation of the targets, comparing to the back-propagation method. 

5 Sparse approximation problem 

The MUSIC algorithm provide very good resolution and separation of targets, 
however it cannot be used to quantify the properties of the targets, such as the 
magnitude of their scattering coefficients. Here we show that, with a little more 
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computational effort, we can also determine the scattering coefficients using a 
sparse image reconstruction approach. 

Let us assume that the imaging domain is discretized as a grid of L D voxels 
(pixels), and ri, I = 1, L D , gives the position of each voxel. Also, we assume 
that pi(u>), 1 = 1, ■■■,L D ', is the scattering factor associated with each voxel in 
the imaging domain. The goal is to find a matrix 

L D 

H{uj) = M^)9(ri^)9{ri,uj) T , (30) 
i=i 

which best approximates the response matrix H(co): 

H(uj) ~ H(lu). (31) 
The above equation can be rewritten as: 

~ 0(w), (32) 
where p(co) is the unknown i^-dimensional vector: 

p(w) = [p 1 {u),...,p L n(u)} T . (33) 
The iV 2 -dimensional vector Q(lu) is obtained by stacking the columns 

H n {uj) = [ H ln (u) ... H Nn (u>) ] T 

of the response matrix H(w): 

e(w) = [ H n (u) ... H ni (uj) ... H 1N (u>) ... H nn (uj) f . (34) 

Also, $(w) is a matrix with N 2 rows and L D columns. Each column 

1 = 1,..., L D , is obtained in a similar way, by stacking the columns of the N x N 

matrix g(r u uj)g(r h uj) T . 

The above system of equations is underdetermined, since the number of 
scattering targets M is much smaller than L D . The common approach to find 
a solution is to consider the equivalent ^-optimization problem [11]: 

p{uj) = argmin ||p(w)|| 2 s.t. ${uj)p{uj) = 6(w), (35) 



where ||p(w)|| 2 = y Sz=i \Pii L0 )\ ■> IS *he Euclidean norm. In this case, the 
unique solution which minimizes the ^-norm, is given by: 

p{w) = <S>Hoj)Q(oj), (36) 

where 

- [S'M^w)]- 1 **^), (37) 
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is the More-Penrose pseudo-inverse of the matrix <J>(w). However, in this case all 
the coefficients of the solution p(u>) are non-zero, and therefore, this solution is 
not in agreement with the fact that the imaging region is sparse, i.e. the number 
of scattering targets is M <C L D . Therefore, this is not the correct solution of 
the approximation problem. In fact, we would like to find the minimum number 
of columns <&;(w) of $(w) which approximate the data vector 6(w). This is a 
sparse approximation problem, and as a measure of sparsity we consider the lo 
norm of p(u>), ||p(w)|| , which simply counts the number of nonzero coefficients 
in the vector p(u>). These non-zero coefficients will give the position, and their 
magnitude will reflect the value, of the scattering coefficients of the targets in 
the L D imaging region. Thus, the sparsest representation requires the solution 
of the Z — optimization problem [12]: 



Unfortunately, this combinatorial optimization problem is NP-hard to solve, 
requiring the enumeration of all possible collections of columns in $(w) and 
searching for the smallest collection which best approximates the data vector 
Q(w). An alternative is the convexification of the objective function, which is 

obtained by replacing the Iq norm with the l\ norm: Hpfw))^ = \Pi( UJ )\- 
The resulting /i-optimization problem: 



is known as Basis Pursuit (BP), and it can be solved using linear program- 
ming techniques whose computational complexities are polynomial [12]. The 
BP method recasts the ^-problem as a linear program, and it has been shown 
that because of the nondiffcrentiability of the li norm, this optimization problem 
leads to unique sparse solutions. However, the BP approach requires the solu- 
tion of a very large convex, nonquadratic optimization problem, and therefore 
still suffers from high computational complexity. For example in a three dimen- 
sional problem, D = 3, with N = 100 transceivers and a discretization grid with 
L = 100, the resulted dimensionality of the <j>(u;) matrix is: 10 4 x 10 6 , and the 
number of unknowns in the vector p(ui) is 10 6 . As an alternative, here we con- 
sider a heuristic approach based on iterative greedy algorithms, which also have 
been proven to give good approximative solutions to the sparse reconstruction 
problem. 

6 Greedy algorithm for sparse approximation 

Matching Pursuit (MP) is a general procedure to compute adaptive signal rep- 
resentations and to extract the signal structure in a given time-frequency dic- 
tionary [13]. Also, it has been shown that the MP algorithm can be used to 
obtain (approximative) sparse solutions of the ^-optimization problem [14-16]. 



p(u>) = argmin ||/o(w)|| s.t. &(uj)p(uj) = Q(ui). 



(38) 



p{to) = argmin H/c^w)!^ s.t. $>(u))p(ui) = Q(w) 



(39) 
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Although the MP algorithm is non-linear, it maintains an energy conservation 
which guarantees its convergence. 

In the case of computational time-reversal imaging, the elements of the time- 
frequency dictionary are given by the columns of the matrix $(w). Using 
this dictionary, the data vector 0(w) can be represented as: 

L D 

0(w) = 5^p,(w)*iM. (40) 
1=1 

The vector 0(w) can be decomposed into: 

e(w) = [$,*(w)*(w)] ||$iMH" 2 $iM + (41) 

where ^(w) is the residual vector after approximating 0(w) in direction of 
Since and ^(w) are orthogonal, we have: 

||*H|| 2 = II9MII 2 - \*Uu)Q(u)\ 2 II^Mir 2 , (42) 

and in order to minimize wc must choose the column &i(u>), such that 

|$;*(o;)0(a;)| ||$j(w)|| 1 is maximum. Thus, starting from an initial approxi- 
mation p(ui) = and a residual %?(cj) = 0(w), the algorithm uses an iterative 
greedy strategy to pick the columns &i(u>) of $>(u>) that are the most strongly 
correlated with the residual. Then, successively their contribution is subtracted 
from the residual, which this way can be made arbitrarily small. The pseudo- 
code of the MP algorithm is: 

1. Initialize the variables: 

T,i<-l,p(w)<-0,#(w)<-e(w). (43) 

2. Find I such that: 

I = arg^max^ |$f H^MIf 1 ■ (44) 

3. Update the estimate of the corresponding coefficient, the residual, and the 
iteration counter: 

- «H + [$f(w)*(w)] ||#iMII" 2 , (45) 

«- - ||$iHir 2 (46) 

i -f— t + 1. (47) 

4. If ||*(w)|| 2 < e ||0(w)|| 2 or f > T then terminate and return /3(a)). Other- 
wise go to 2. 
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The stopping criterion in the step 4 requires the residual to be smaller than 
some fraction < e <C 1 of the data vector. Also, the computation stops if 
the number of iterations exceed the maximum number allowed T. Although 
the asymptotic convergence of MP algorithm can be easily proven, the resulting 
approximation after any finite number of steps will in general be suboptimal. In 
the case of noise, the MP algorithm is used to obtain an approximative sparse 
solution by simply stopping the iteration when the projection of the residual on 
the chosen direction becomes smaller than a threshold < r < 1 [14-16]. 

Thus, in the case of noise, the MP algorithm is modified simply by replacing 
the stopping condition with: |<J>;*(u;)\I>(u;)| 2 ||$z(w)||~ 2 ||*(o))|p 2 < r. After the 
computation is finished, the image is formed by plotting pi(u>) at location r/, 
I = 1,...,L D . We should note that the algorithm works in both cases, when 
pi(u>) is complex, or when only its magnitude |jO;(w)| is given. In the later case, 
one should plot the absolute values of the computed coefficients \pi{to)\. 

7 Implementation and numerical results 

It is important to note that one can implement the MP algorithm such that the 
elements of the matrix $ do not need to be stored. In fact, one can compute 
the columns of <E> at every step of the algorithm, since the Green function is 
known. Thus, the algorithm requires only operations with vectors of length N 2 , 
which is feasible on standard personal computers. However, since the number of 
vector- vector multiplications per iteration step is high, L 2 , a parallel implemen- 
tation of the algorithm is desirable. We have implemented the MP algorithm 
for the NVIDIA GPU (Graphics Processing Unit) platform. Recently, NVIDIA 
has released a general purpose oriented API for its graphics hardware, called 
CUDA [17]. In addition, NVIDIA has developed CUBLAS which is a GPU 
optimized version of BLAS library (Basic Linear Algebra Subroutines) built on 
top of CUDA [18]. The newly developed GPUs now include fully programmable 
processing units that follow a stream programming model and support vector- 
ized single and double precision floating-point operations. For example, the 
CUDA computing environment provides a standard C like language interface 
to the NVIDIA GPUs. The computation is distributed into sequential grids, 
which are organized as a set of thread blocks. The thread blocks are batches 
of threads that execute together, sharing local memories and synchronizing at 
specified barriers. CUBLAS library provides functions for: (i) creating and de- 
stroying matrix and vector objects in GPU memory; (ii) transferring data from 
CPU mainmcmory to GPU memory; (iii) executing BLAS on the GPU; (iv) 
transferring data from GPU memory back to the CPU mainmemory. CUBLAS 
defines a set of fundamental operations on vectors and matrices which can be 
used to create optimized higher-level linear algebra functionality: (i) Level 1 
BLAS perform scalar, vector and vector-vector operations; (ii) Level 2 BLAS 
perform matrix-vector operations; (iii) Level 3 BLAS perform matrix-matrix 
operations. However, in its current version CUBLAS does not offer direct sup- 
port for operations involving vectors and matrices of complex numbers. In the 
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case of the MP algorithm, one can easily overcome this drawback by storing 
separately the real and imaginary part of the vectors and matrices. This way, 
complex vector-vector and matrix-vector operations can be reduced to oper- 
ations involving only real numbers. The parallel implementation of the MP 
algorithm requires only Level 1 or Level 2 (if <!> is stored) BLAS operations. 
Our numerical tests have shown that the parallel GPU implementation versus 
a standard CPU BLAS implementation of the MP algorithm reaches a speed 
up of a maximum of 31 times in single precision and respectively 21 times in 
double precision [19]. 

Let us now consider a two dimensional scenario, consisting of N = 100 
linearly distributed transceivers, separated by d = A/2 and located at x n = 
[ n\/2 + a/2 — NX/ '4 ] , where a = 100A is the side of the imaging area. 
The computational image grid was also set to L x L = 100 x 100. This scenario 
will generate a dictionary $ of size N 2 x L 2 = 10 4 x 10 4 and an unknown 
vector p of size A 2 = 10 4 . The number of targets is set to M = 5 and their 
position is randomly generated in the imaging area. The scattering coefficients 
are generated randomly from a uniform distribution such that 1 < p m < 10, 
m = 1, ...,M. We consider both cases, with and without noise. The signal to 
noise level is set to SNR — 2, as for the described MUSIC algorithm case. 

First, let us discuss the case without noise. The threshold parameter and the 
maximum number of iterations were set to s = 10~ 4 , and respectively T — N. 
In Figure 4 we give the initial and the computed arrangement of the targets, 
and their initial and computed values of the scattering coefficients. One can 
see that the agreement between the initial values and the computed ones is 
very good, and the MP algorithm can solve the problem almost perfectly. In 
Figure 5 we give the real and imaginary part of the initial and reconstructed 
response matrix. The error is less than 1% in both cases. Also, Figure 6 shows 
the computed image using the MP algorithm. One can see that the peaks are 
very sharp and their position and amplitude reflects correctly the initial position 
and the magnitude of the scattering coefficients. This example shows that the 
MP algorithm can solve almost perfectly (in the limits of the image resolution) 
the computational time-reversal imaging problem if the response matrix is not 
affected by noise. 

Let us consider the case with noise, using the same arrangement of the 
targets, with the same values for the scattering coefficients. In this case the 
threshold r plays a very important role in the MP algorithm. We consider two 
different values r = 10~ 3 and r = 10~ 2 . In Figure 7 we give the real and the 
imaginary part of the initial and reconstructed response matrix. One can see 
that the reconstructed response matrix contains a much lower amount of noise 
than the original response matrix. Also, by increasing the value of r, the amount 
of noise in the reconstructed response matrix decreases dramatically. This can 
be seen also on the computed images, given in Figure 8. The amount of noise in 
the computed images is very small, given the SNR = 2 level in the perturbed 
response matrix. Also, the position of the peaks, corresponding to the targets, 
and their magnitude are still very well maintained. This means that the sparse 
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approximation, given by the MP algorithm, and the threshold parameter r , 
acts as a denoising method, and respectively as a denoising parameter, and it 
can be succesfully used in computational time-reversal imaging. 

8 Conclusion 

We have presented a sparse approximation method for computational time- 
reversal imaging. The method is formulated entirely in the frequency domain, 
and it is based on an adapted version of the matching pursuit algorithm, which 
can be successfully used to compute an accurate sparse approximation of the 
frequency response matrix. This approach can be used for denoising the com- 
puted time-reversal images, and to determine the magnitude of the scattering 
coefficients of the targets embedded in homogeneous media, in the presence of 
moderate to high noise levels. Also, in comparison to the back-propagation and 
the null subspace projection methods, the described approach provides a better 
resolution. However, the sparse approximation method is computationally more 
expensive than the traditional approaches. 
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Figure 1: Geometry for the time-reversal imaging experiment, containing N 
transceivers and M scattering targets. 
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Figure 4: Sparse approximation of time-reversal image in case of M = 5 targets: 
(a) the initial location of the targets and their scattering coefficients; (b) the 
computed location of the targets and their scattering coefficients. 
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(a) 
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Figure 5: Sparse approximation of time-reversal image in case of M — 5 targets 
without noise: (a-b) the real and imaginary part of the initial response matrix 
(c-d) the real and imaginary part of the computed response matrix. 
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(a) (b) 

Figure 6: The computed time- reversal image, using the sparse approximation 
approach, in case of M = 5 targets, without noise. 
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(e) (f) 



Figure 7: Sparse approximation of time-reversal image in case of M = 5 targets, 
with noise: (a-b) the real and imaginary part of the initial response matrix with 
the noise level SNR = 2; (c-d) the real and imaginary part of the computed re- 
sponse matrix with r = 1CP 3 ; (e-f) the real and imaginary part of the computed 
response matrix with r = 1CP 2 . 
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